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We find diffraction-free beams for graphene and MoS2-type honeycomb optical lattices. The 
resulting composite solutions have the form of multi-vortices, with spinor topological charges (n, 
nil). Exact solutions for the spinor components are obtained in the Dirac limit. The effects of 
the valley degree of freedom and the mass are analyzed. Passing through the van-Hove singularity 
the topological structure of the solutions is modified. Exactly at the singularity the diffraction-free 
beams take the form of strongly localized one-dimensional stripes. 


The dynamics of an optical wave can be engineered 
by either modulating the shape of the wave itself or by 
changing the properties of the medium. Over the past 
years a large amount of effort has been devoted in both 
of these directions. Of central importance is the ability 
to control the inherent diffraction of propagating waves. 
One limiting case is to completely suppress diffraction 
leading to the generation of the so-called diffraction- 
free beams (DFBs). In homogeneous media there are 
two main classes of DFBs: The first one is the Bessel 
beam mm which, in terms of ray optics, can be de¬ 
scribed as an extended and elongated focus. The second 
class is the Airy beam mm which, in terms of ray optics, 
can be described as a caustic. Along these two directions 
a variety of different classes of beams have been proposed 
and observed mm. a diffraction-resisting hybrid having 
the form of the Bessel beam but being able to trans¬ 
versely accelerate along pre-defined trajectories has also 
been realized muni. 

The other ingredient which can be appropriately engi¬ 
neered is the medium. Periodic variations of the refrac¬ 
tive index can be easily created, using for example opti¬ 
cal induction m, resulting to a host of discrete wave 
phenomena [12]. DFBs can be constructed in square 
lattices by utilizing the Floquet-Bloch structure of the 
system m • The equivalence between the equations gov¬ 
erning the light dynamics in these arrays and the ones in 
solid state physics or in quantum theories, offers the op¬ 
portunity to explore in optical settings phenomena pre¬ 
dicted in branches of science beyond optics. Perhaps, 
the most eminent representative of such a structure is 
graphene, which was recently made available p3]. Effort 
has been made to design graphene-like materials with 
modified properties. In M 0 S 2 m a minigap opens at 
the Dirac points, a property highly desirable for poten¬ 
tial electronic applications. A bandgap can also open 
up by deforming a honeycomb lattice p~6HT8] . Relevant 
phenomena have been studied in honeycomb waveguide 
arrays including conical diffraction m , the trembling 
motion (Zitterbewegung) of electrons [20 , the Klein tun¬ 
neling m , topological insulators m and the relation 
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between pseudospin and orbital angular momentum J21 . 

In this work, we find DFBs for graphene and M 0 S 2 - 
type honeycomb optical lattices. We introduce a physi¬ 
cally relevant spinor field decomposition method to ana¬ 
lyze the underlying structure of these waves. The result¬ 
ing composite solutions have the form of multi-vortices 
or semi-vortices, in the sense that the two spinor compo¬ 
nents are associated with different topological charges (n, 
n± 1). Asymptotically, close to the Dirac point, the sys¬ 
tem is described by the Dirac-Weyl equation (massless or 
massive). In this limit, we analytically obtain diffraction- 
free composite multi-vortex solutions which are in agree¬ 
ment with our numerical results. Furthermore, the effect 
of the valley inequivalence of the K and K' points leads 
to vorticies with different topological charges. As the 
propagation constant traverses the band structure (BS), 
the diffraction-free solutions undergo a transition exactly 
at the van-Hove singularity. As a result, after this transi¬ 
tion, the spinor components of the DFBs have the same 
vorticity. Exactly at the van-Hove singularity three dif¬ 
ferent solutions are obtained which have the form of an 
array of infinite extent in one direction but strongly lo¬ 
calized in the orthogonal direction. We also discuss ex¬ 
tensions of this work to generic lattices with more than 
one “atoms” per unit cell. 

We consider the paraxial dynamics of light waves in 
normalized coordinates iip z + (1/2)+ V(r) 2 p = 0, 
where ^ is the field amplitude, r = (x,y) are the trans¬ 
verse coordinates and z is the propagation coordinate, 
vi = d* + di is the diffraction operator and V(r) is 
the potential that is proportional to the refractive index 
contrast. In particular, we select a honeycomb periodic 
potential of the form V(r) = (Vo/c§) | 22 j =1 — 

c 0 | 2 , where G x = Gx, G = 2tt/L, Gj = W~ 1 (ti/3)G 1 , 
L is the lattice period, and R the rotation matrix and 
<f> 2 j+i = x/2, <f> 2 j = — x/2. In the rest of the paper 
we choose Co = 6 and % = 0 for graphene lattices and 
X = 7r/20 for optical M 0 S 2 lattices (honeycomb with 
sublattice symmetry breaking) [see Fig. [I]. Note that 
similar lattices can be generated via optical induction as 
shown in [22]. Due to symmetry breaking, in M 0 S 2 lat¬ 
tices a bandgap opens up at the K and K' points of the 
Brillouin zone [Fig. 00 )]- 
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FIG. 1. Typical index profiles of a graphene (a) and a M 0 S 2 
(b) lattice (darker areas represent higher index); (c) lattice 
structure; (d) contour plot of the graphene BS, (blue curve is 
the van-Hove singularity); (e) BS for graphene (dashed red) 
and M 0 S 2 (solid black) lattices. 



FIG. 2. A diffract ion-free composite semi-vortex with topo¬ 
logical charges Q = (0, —1) supported by a graphene type 
lattice for Vo = 15 and L = 1.81. The propagation constant 
is /3 = —19.099 and the Dirac point is fo = —18.905. In 
the first and second rows the intensity and phase patterns 
are shown. In column (a) the total field ^ is shown, while in 
columns (b) and (c) the isolated spinor components and 
b are presented (the insets show the corresponding sublat¬ 
tices). The inset of (a) is a typical Bloch mode amplitude 
along the trajectory shown in inset (d). 

For relatively high index constrasts, we can apply 
coupled mode theory (CMT). We denote the lowest 
eigenmode of an isolated waveguide of the potential by 
U?n,n (*G y) 1 where the subscripts (ra, n) label the trans¬ 
lations over the unit cell [Fig. flTc)] , while the super¬ 
script j = {A, B} = {1,2} denotes the two sites in 
each unit cell. We decompose the optical field as ip = 
'Em,n,j 4n,n( z ) u L,n( r )i where c? m n are the ^-varying 
amplitudes which can be written in spinor form as c m ^ n = 
(c^ n ,c^ n ). Following the relevant calculations we ob- 

tain i£4n, n + t T,(m,n) C m,l + (-l) i+1 /?0 4n,n = 0 where 
t is the coupling strength between adjacent waveguides, 
(m, n) denotes coupling over first neighbors, and ±fo are 
the relative shifts of the propagation constants due to the 
sublattice asymmetry resulting to a bandgap width 2fo. 

The continuous limit of the CMT equations can be 
derived [c m ^ n (z) —>■ u(x,y,z) = (a, b)] by selectively ex¬ 
citing the lattice with a broad beam in the vicinity of par¬ 
ticular k -vectors of the BS [23]. Specifically, near the ver¬ 
tices (K and K r ) of the 1st Brillouin zone of Fig. OF) the 
following Dirac equation is derived id z a — D(dfo)-\-foa = 



FIG. 3. Same as in Fig. 1 for a diffraction-free composite 
multi-vortex with topological charges Q = (-1,-2). 

0 id z b + D(d±a) — fob = 0, where D = y/3tL/ 2 and 
fo describes the effective mass in the context of M 0 S 2 
and is zero for graphene. In addition d± = d x ± idy = 
e ±uf> ($ r ± and the signs refer to the two inequiva¬ 

lent Dirac points (upper sign is K and lower sign is K r ). 
Without loss of generality, we set D = 1 for the rest of 
the paper. 

We have found DFB solutions of the Dirac system 
which are expressed in terms of Bessel functions as 

a(r,<f>,z) = J n {\r)e in<t ’e~ i ^ z , (1) 

b(r, </>, z) = gJ n ±i(Xr)e^ n±1 ^e-^ z , (2) 

where g(/3) = T\W+M/W Zr Po), and A = - 0$. 

Thus the spinor u = (a, b) consists of two vortices with 
unequal topological charges Q =» (n, n =b 1). We call this 
DFB a composite multi-vortex DFB. We note that for 
graphene lattices fo = 0 and thus the two spinor com¬ 
ponents are equally “inhabited” (\g\ = 1). On the other 
hand, the presence of an effective mass fo gives rise to a 
completely different behavior. As we approach the band 
edges /3 —>> fo (/3 — fo) the relative amplitude of b (a) 
goes to zero and we end up with a single vortex with vor- 
ticity n (nil), respectively. As we move away from the 
gap edges, the relative amplitude of the two components 
become equal. 

It is interesting to point out that, in contrast to the 
Dirac equations, within the paraxial framework both 
spinor components are described by a single wavefunc- 
tion [2T] . Thus, it is physically relevant to follow a similar 
decomposition, for the paraxial model. In particular we 
separate the wavefunction of the paraxial equation dis¬ 
cretely according to its spatial location as 'ipA(Rm,n^ z ) 
and z )- I n above equations R J m , n are the 

position vectors of the lattice elements A and B [see 
Fig. 0 c)] located in the same Wigner-Seitz cell (m,n). 
Consequently, there is absolutely no meaning in any 
form of interference between these two components of the 
Dirac equation (at least when the CMT is applicable). 

The paraxial equation supports Floquet-Bloch modes 
having the form u^fo) = Vk{r)e Lk ‘ r ~ l ^^ z where Vk(r) 
is a function with the same periodicity as the lattice. 
However, if u & is a Floquet mode the same holds for 
for an arbitrary phase £. In the case of a single Floquet 
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FIG. 4. Valley and mass effects of M 0 S 2 composite semi- 
vorticies for /3 = —20.424, Vo = 15, and L = 1.81 encircling 
the K ( K') point in second (third) row leading to Q = (—1, 0) 
[Q = (1,0)], respectively. The intensity pattern shown in the 
first row is identical in both cases. 

mode, this might seem a trivial generalization, however, 
when superimposing Floquet modes their relevant phase 
becomes important. We make two different choices for 
this phase factor and denote the corresponding Floquet 
modes as u J k , j = {A, B}. In particular, we choose all the 
modes to have a positive real amplitude on the sublattice 
with (ra,n) = (0,0) (u k (R^ 0 ) > 0 and u k (R B 0 ) > 0). 

The function /3(k) describes the BS surface. As we 
can see in Fig. [ljd) the isocontours (/3 = constant) sur¬ 
round the K and K' or the T points of the BS and, as 
a result, can be described in terms of a polar coordinate 
ip. The separating curve corresponds to the van-Hove 
singularity. We can invert the relation /3(k) = /3q and 
solve in terms of the Bloch momentum as k(cp,/3o) for 
each curve. We conclude that the general form of a DFB 
is = / 0 2 " A( K (p)e lG ^u^ k ^e' lk ^' r ~ l/3 ^ z d(p for arbitrary 
values of the amplitude A(ip) and the phase 0((p). Note, 
that, as we described, for the same value of {3 surround¬ 
ing the Dirac cones there are four different families of 
diffraction free beams depending whether we encircle K 
or K' and on the Floquet modes u k or u k . 

In Fig. [2] we show a DFB that is generated by equal 
(A(ip) = 1) in-phase (0(<p) = 1) excitation of the Flo¬ 
quet modes u k along the, almost circular, 1st band closed 
contour that surrounds the K' points. Thus, we expect 
that the asymptotic solutions [Eqs. 0-0] are in good 
agreement with the numerical solutions. The underly¬ 
ing structure of the wave is revealed by decomposing the 
field to its spinor components as we previously described. 
Numerically, in each unit cell (m, n) we calculate the op¬ 
tical field at the index maxima of the sites A and B. 
We then set the value of the field everywhere inside this 
cell to be equal to this latter calculated value. Since the 
Floquet modes around the Dirac cone consist of a pe¬ 
riodic array of vortices, it is relevant to keep only one 



FIG. 5. Directional 1st band DFB supported by the van-Hove 
singularity when traversing (a) the closed hexagonal path and 
(b) two parallel line segments. The structure of the modes is 
shown in the right insets of (b). A graphene lattice is chosen 
with Vo — 15 and L = 1.81, while [3 — —19.442. 



FIG. 6. DFB encircling the T-point for (a) /3 — —19.468 and 
(b) P — —19.784 in a graphene lattice with Vo = 15 and 
L = 1.81. 

out of three unit cells to avoid the imaging of this local 
vortex structure. Using this decomposition, we clearly 
see the formation of a composite semi-vortex with topo¬ 
logical charges Q = (0,-1) in agreement with Eqs. 0 - 
© with n = 0 for the K'. In Fig. [ 3 ] we choose the 
same parameters as in Fig. [2] but we impose a vortex 
phase structure O(cp) = —27 d(ip)/L, where l(ip) is the 
arclength. This leads to a composite multi-vortex DFB 
with Q = (—1, —2) [see Eqs. 0 -© with n = —1]. 

In Fig. [4] we analyze the effects of (i) the valley de¬ 
gree of freedom and (ii) the mass on the formation of 
DFBs. We select 1st band contours with the same (3 
encircling the inequivalent K and K' points. A typi¬ 
cal Floquet-mode beam intensity is localized on sublat¬ 
tice B [Fig. EJa) inset] although the index is higher on 
sublattice A [Fig. [ljb)] . This happens due to the local 
vortex structure that enhances the energy of the respec¬ 
tive Hamiltonian. Thus, we choose the Floquet modes 
u k which constructively interfere on sublattice B and 
A(ip) = 1, O(ip) = 0. Due to the valley degree of free¬ 
dom (inequivalence of K and K') the generated compos¬ 
ite semi-vortices, although they have the same intensity 
pattern, are associated with charges Q = (—1,0) and 
Q = (1,0), as also predicted from Eqs. ([l])-([2|. The effect 
of the mass leads to another interesting feature. Specifi¬ 
cally, the resulting composite DFB has strongly unequal 
amplitudes in its spinor components as we approach the 
band edges. In Fig. [4] the relative maximum amplitude 
of the spinor components is max \ ^ A \/ max \ 'ip B \ = 0.06. 

As we move further away from the Dirac points the so¬ 
lutions remain topologically equivalent to those presented 
above, having the same topological charges Q. Eventu¬ 
ally, we reach the van-Hove singularity which is the con¬ 
tour that passes through the boundary midpoints of the 
hexagonal 1st Brillouin zone. In Fig. [5j[a) we equally ex¬ 
cite the path that constitutes the van-Hove singularity. 
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As we can see the DFB is associated with three dominant 
directions. The structure of the DFBs at the singularity 
is revealed by selectively exciting only two parallel lines 
segments. In Fig. [5^b) we clearly see a DFB that is of 
infinite extent in the ^/-direction and strongly confined 
along the x-direction. Three such solutions exist with 
a 27r/3 rotational symmetry. Within the framework of 
the CMT these solutions are expressed analytically as 
Cl-I = (—1)*k(1, —(A) + /3±)/t)e~ ll3z and c m , n = 0 oth¬ 
erwise, where /3± = ± \/3^ + t 2 [the indices are shown 
in Fig. file)]. The /3_ (/3 + ) eigenvalue lies in the first 
(second) zone of the BS [see top (bottom) right inset of 
Fig. 5jb)]. In the graphene limit (/3q = 0) the nonzero 
amplitudes of the modes in the two sublattices become 
equal. In M 0 S 2 lattices when fo/t 1 the amplitudes 
in one sublattice become negligible resulting to a single 
out-of-phase straight array of diffraction-free light beams. 

The van-Hove signularity separates the BS into topo¬ 
logically inequivalent regions. Passing through this sin¬ 
gularity we approach the neighborhood of the T point 
of the Brillouin zone (which the DFB encircles) and the 
Floquet modes lose the local vortex structure. Typical 
DFBs in this case are shown in Fig. [6j Specifically, in 
Fig. [6](a) the path is close to the van-Hove singularity 
and exhibits a “starfish” like structure which is reminis¬ 
cent of the directional modes shown in Fig. |5j Moving 


towards the T point the integration path as well as the 
DFB take a cicular form [Fig. [6](b)] . 

The analysis presented above is not limited to honey¬ 
comb lattices but can be applied to any lattice with more 
than one “atom” per unit cell. Our investigations show 
that the presence of the van-Hove singularity separates 
regions where the composite DFBs have different topo¬ 
logical charges. However, in many lattices the van-Hove 
singularity exists exactly at the band edge. For exam¬ 
ple, in the case of a square diatomic lattice (where the 
singularity exists in the band edge) both spinor compo¬ 
nents have the same vorticity everywhere inside the band. 
Also, exactly at the singularity elongated DFBs always 
exist. 

Finally, we would like to point out that DFBs not only 
constitute a fundamental class of solutions and the build¬ 
ing blocks for other classes of beams but are also impor¬ 
tant in unveiling the fundamental properties of the sys¬ 
tem. For example, conical diffraction can be explained 
by expanding the initial condition to the diffraction-free 
beams supported by the system [2Tj . 
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